proc ds2;
  package math / overwrite = yes;
    method radians(double angle) returns double;
      return angle * constant('PI') / 180;
    end;

    method rational(double x, double maxden, in_out integer p, in_out integer q) returns double;
      /*
      ** FROM: https://www.ics.uci.edu/~eppstein/numth/frap.c
      ** FROM: https://stackoverflow.com/questions/95727/how-to-convert-floats-to-human-readable-fractions
      **
      ** find rational approximation to given real number
      ** David Eppstein / UC Irvine / 8 Aug 1993
      **
      ** With corrections from Arno Formella, May 2008
      **
      ** Modified for Proc DS2, Richard DeVenezia, Jan 2020.
      **
      ** usage: rational(r,d,p,q)
      **   x is real number to approx
      **   maxden is the maximum denominator allowed
      **   p is return for numerator
      **   q is return for denominator
      **   returns 0 if no problems
      **
      ** based on the theory of continued fractions
      ** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
      ** then best approximation is found by truncating this series
      ** (with some adjustments in the last term).
      **
      ** Note the fraction can be recovered as the first column of the matrix
      **  ( a1 1 ) ( a2 1 ) ( a3 1 ) ...
      **  ( 1  0 ) ( 1  0 ) ( 1  0 )
      ** Instead of keeping the sequence of continued fraction terms,
      ** we just keep the last partial product of these matrices.
      */

      declare integer m[0:1,0:1];
      declare double startx e1 e2;
      declare integer ai t result p1 q1 p2 q2;

      startx = x;

      /* initialize matrix */
      m[0,0] = 1; m[1,1] = 1;
      m[0,1] = 0; m[1,0] = 0;

      /* loop finding terms until denom gets too big */
      do while (1);

        ai = x;

        if not ( m[1,0] * ai + m[1,1] < maxden ) then leave;

        t = m[0,0] * ai + m[0,1];
        m[0,1] = m[0,0];
        m[0,0] = t;

        t = m[1,0] * ai + m[1,1];
        m[1,1] = m[1,0];
        m[1,0] = t;

        if x = ai then leave;     %* AF: division by zero;

        x = 1 / (x - ai);

        if x > 2147483647 /*x'7FFFFFFF'*/ then leave;  %* AF: representation failure;
      end;

      /* now remaining x is between 0 and 1/ai */
      /* approx as either 0 or 1/m where m is max that will fit in maxden */
      /* first try zero */

      p1 = m[0,0];
      q1 = m[1,0];
      e1 = startx - 1.0 * p1 / q1;

      /* now try other possibility */

      ai = (maxden - m[1,1]) / m[1,0];

      m[0,0] = m[0,0] * ai + m[0,1];
      m[1,0] = m[1,0] * ai + m[1,1];

      p2 = m[0,0];
      q2 = m[1,0];
      e2 = startx - 1.0 * p2 / q2;

      if abs(e1) <= abs(e2) then do;
        p = p1;
        q = q1;
      end;
      else do;
        p = p2;
        q = q2;
      end;

      return 0;
    end;
  endpackage;
run;

  package pen / overwrite = yes;
    declare private int down;
    declare private varchar(20) pencolor;
    declare private varchar(20) fillcolor;
    declare private int size;

    method up();
      down = 0;
    end;
    method down();
      down = 1;
    end;
    method isDown() returns int;
      return down;
    end;
    method isUp() returns int;
      return not down;
    end;
    method color() returns varchar(20);
      return pencolor;
    end;
    method color(char(8) color);
      pencolor = color;
    end;
    method size() returns int;
      return size;
    end;
    method size(int newsize);
      size = newsize;
    end;
  endpackage;
run;

  package canvas / overwrite=yes;
    declare private int xpixels  ypixels;
    declare private varchar(41) statements_table;
    declare private varchar(300) path;
    declare private varchar(300) filename;
    declare private varchar(41) annoname;
    declare private package pen pen;

    declare package sqlstmt q;
    declare varchar(1000) s;

    method canvas(
      varchar(41) statements_table,
      varchar(300) path,
      varchar(300) filename
    );
      canvas(1200, 675, statements_table, path, filename); * 16:9;
    end;
    method canvas(
      int xpixels,
      int ypixels,
      varchar(41) statements_table,
      varchar(300) path,
      varchar(300) filename
    );
      declare varchar(41) annoname;
      annoname = translate(filename,'__',' .');

      canvas (xpixels, ypixels, statements_table, annoname, path, filename);
    end;

    method canvas(
      int xpixels,
      int ypixels,
      varchar(41) statements_table,
      varchar(41) annoname,
      varchar(300) path,
      varchar(300) filename
    );
      this.statements_table = statements_table;
      this.annoname = annoname;
      this.path = path;
      this.filename = translate(filename,'_','.');

      this.xpixels = xpixels;
      this.ypixels = ypixels;

      q = _new_ sqlstmt('insert into ' || statements_table || ' values (?)', [s]);

      s = 'data ' || annoname || '; %sgvardcl; keep drawspace;';
      q.execute();
      s = 'drawspace=''graphpixel'';';
      q.execute();
    end;
    method translate(in_out double x, in_out double y);
      %* map turtle canvas x,y to pixel image x,y;
      %* 0,0 is center of canvas;
      %* 0,0 is bottom left of pixel image;
      x + xpixels / 2;
      y + ypixels / 2;
    end;
    method rgb_to_color(int r, int g, int b) returns char(8);
      return 'CX' || put(mod(r,256),HEX2.) || put(mod(g,256),HEX2.) || put(mod(b,256),HEX2.);
    end;
    method setpen(package pen pen);
      this.pen = pen;
    end;
    method line(double x1, double y1, double x2, double y2);
      if pen.isUp() then 
        return;

      translate (x1,y1);
      translate (x2,y2);

      s = '%sgline(x1='||round(x1)||
                 ',y1='||round(y1)||
                 ',x2='||round(x2)||
                 ',y2='||round(y2)||
                 ',linecolor='''||pen.color()||''''||
                 ',linethickness='||pen.size()||
                 ');';

      q.execute();
    end;
    method arc(double xc, double yc, double radius, double heading, double extent,int steps);
      declare double x1 y1 x2 y2;
      declare double p1 q1 p2 q2;
      declare double arcstep;

      s = 'note="Arc-start"; keep note;output;';
      q.execute();

      x1 = xc + radius * cos (heading * constant('pi') / 180);
      y1 = yc + radius * sin (heading * constant('pi') / 180);

      p1 = x1;
      q1 = y1;

      translate (p1,q1);

      if radius < 0 then
        extent = -extent;

      if steps = 0 then
        steps = abs(extent / 2);

      arcstep = extent / steps;

      do heading = heading + arcstep to heading + extent by arcstep;
        x2 = xc + radius * cos (heading * constant('pi') / 180);
        y2 = yc + radius * sin (heading * constant('pi') / 180);

        p2 = x2;
        q2 = y2;

        translate (p2, q2);
  
        s = '%sgline(x1='||round(p1)||
                   ',y1='||round(q1)||
                   ',x2='||round(p2)||
                   ',y2='||round(q2)||
                   ',linecolor='''||pen.color()||''''||
                   ',linethickness='||pen.size()||
                   ');';
        q.execute();

        x1 = x2;
        y1 = y2;

        p1 = p2;
        q1 = q2;
      end;

      s = 'note="Arc-end";output;';
      q.execute();
    end;
    method filename() returns varchar(300);
      return filename;
    end;
    method write(double x, double y, varchar(200) text, varchar(400) options);
      declare double w;

      if pen.isUp() then 
        return;

      translate (x,y);

      s = '%sgtext (x1='||round(x)
              ||  ',y1='||round(y)
              ||  ',label='||quote(text)
              ||  ',width='||xpixels
              ||  ',widthunit="pixel"'
              ||  ifc(length(options), ','||options, '')
              ||  ');';
      q.execute();
    end;
    method close();
      s = '%fdelete('||path||'/'||filename||'.png)';
      q.execute();

      s = 'ods listing gpath=' || quote(path) || ';';
      q.execute();

      s = 'ods graphics /' ||
          ' outputfmt=png' ||
          ' imagename="' || filename || '"' ||
          ' reset=index' ||
          ' width=' ||round(xpixels)||'px' ||
          ' height='||round(ypixels)||'px' ||
          ';'
      ;
      q.execute();

      s = 'proc sgrender sganno='||annoname||' template=DS2GTI;run;';
      q.execute();
      q.delete();
    end;
  endpackage;
run;

  package turtle / overwrite=yes;
    declare private double x y heading;
    declare private package canvas canvas;
    declare private package math math();

    declare private package pen pen;

    forward home;

    method turtle(package canvas canvas);
      this.canvas = canvas;
      pen = _new_ pen();
      canvas.setpen(pen);
      pen.up();
      home();
    end;
    method log();
      put 'NOTE:' x= y= heading=;
    end;
    method advance(double distance);
      declare double xnew ynew;

*put 'NOTE: advance ' distance=;

      xnew = x + distance * cos(math.radians(heading));
      ynew = y + distance * sin(math.radians(heading));

      canvas.line(x,y,xnew,ynew);

      x = xnew;
      y = ynew;
    end;
    method retreat(double distance);
      advance (-distance);
    end;
    method right(double angle);
      heading = mod (heading - angle, 360);
    end;
    method left(double angle);
      heading = mod (heading + angle, 360);
    end;
    method setxy(double xnew, double ynew);
      canvas.line(x,y,xnew,ynew);
      x = xnew;
      y = ynew;
    end;
    method setx(double xnew);
      canvas.line(x,y,xnew,y);
      x = xnew;                                                     
    end;
    method sety(double ynew);
      canvas.line(x,y,x,ynew);
      y = ynew;
    end;
    method getx() returns double;
      return x;
    end;
    method gety() returns double;
      return y;
    end;
    method heading(double heading);
      this.heading = heading;
    end;
    method heading() returns double;
      return heading;
    end;
    method home();
      x = 0;
      y = 0;
      heading = 0;
    end;
    method arc(double radius, double extent, int steps);
      %* radius > 0, center of arc is on heading left normal - arc direction is counter clockwise;
      %* radius < 0, center of arc is on heading right normal - arc direction is clockwise;
      %* how much of the arc to draw;
      %* steps, not used;

      %* turtle heading is tangent to arc to be drawn;
      declare double xc yc;

      xc = x + radius * cos(math.radians(heading + 90)); %* +90, left normal;
      yc = y + radius * sin(math.radians(heading + 90));

      canvas.arc(xc,yc,radius,heading-90,extent,steps);

      x = xc + radius * cos(math.radians(heading + extent - 90));
      y = yc + radius * sin(math.radians(heading + extent - 90));

      heading = mod (heading + extent, 360);
    end;
    method dot();
    end;
    method down();
      pen.down();
    end;
    method up();
      pen.up();
    end;
    method color(char(8) color);
      pen.color(color);
    end;
    method size(double size);
      pen.size(size);
    end;
    method size() returns double;
      return pen.size();
    end;
    method write(varchar(200) text, varchar(400) options);
      canvas.write(x,y,text,options);
    end;
  endpackage;
run;
quit;